' PRINT "Transducer directivity and pdf    3/2/91 rv 4/27/93 csc "

DIM wr(8000),b(8000),bm(100),wdb(100),phi(100)

'use transformation of pdf to compute w(b) for x = ka sin(t)
'D(x) = 2J1(x)/x
'reduced wr(t) (ka)^2 = (sin(t))^2/(4cos(t) (J0(x) - D(x)) )
' w(t) = (ka)^2 wr(t)

 pi = 4*ATN(1)
 lgcv =20/LOG(10)
 lg10 = LOG(10)
 da = lg10/10 :'10 steps per decade
 mmax = 12
bmin = .07
 esp = .0000001 :'set level for the dB calc and glitches
  upL = 5E+09 :'set max amplitude parameter in Bessel calc
  
   'coef for J0(x)

 jc2=-2.2499997#: jc4=1.2656208#: jc6=-.3163866: jc8=.0444479: jc10=-.0039444
 jc12 =.00021
 jd0=.79788456#: jd1=-.00000077#: jd2=-.0055274: jd3=-.00009512#
 jd4=.00137237#: jd5=-7.2805E-04: jd6=.00014476#
 je0=-.78539816#:je1=-.04166397#:je2=-.00003954#:je3=.00262573#
 je4=-.00054125#:je5=-.00029333#:je6=.00013558#
  
   'coef for J1(x)
 
 ja0=.5: ja2=-.56249985#: ja4=.21093573#: ja6=-.03954289#: ja8=.00443319#
 ja10=-.00031761#: ja12=.00001109#
 jb0=.79788456#:jb1=.00000156#: jb2=.01659667#: jb3=.00017105#: jb4=-.00249511#
 jb5=.00113653#:jb6=-.0002033
 jg0=-2.3561945#: jg1=.12499612#:jg2=.0000565: jg3=-.0063788:jg4=.00074348#
 jg5=.00079824#: jg6=-.00029166#
 

 cyl$ = " circular piston"

 n1 = 1
nmax = 2000
ka2 = 10
tmax = 60
dt = tmax/nmax

 PRINT" Directional response for circular piiston transducer"
 PRINT"Radius ka"

11   'menu
     
     PRINT"input 'f' for ka"
     PRINT"      't'  number of steps in theta"
     PRINT"      'pp' print parameters"
     PRINT"input 'c' to compute"
     
     PRINT"      'mf' to make file"
     PRINT"      'm' for menu"
     PRINT"      'q' to quit"
12  PRINT "go  'm' for menu";:INPUT q$
    
    
     IF q$ = "f" GOTO 60
     IF q$ = "t" GOTO 90 

     IF q$ = "pp" GOTO 100
     IF q$ = "c" GOTO 200

     IF q$ = "mf" GOTO 5000
     IF q$ = "q" GOTO 9000
     IF q$ = "m" GOTO 11
     GOTO 12

60 PRINT"old ka =";ka2;"new =";:INPUT ka2
 GOTO 12

90 PRINT" theta goes 0 t0 ";tmax ;" ,deg  new tmax =";:INPUT tmax
 PRINT "old number of steps =";nmax;" new steps in theta =";:INPUT nmax
 dt = tmax/nmax
PRINT "old n1 =";n1;" new n1 =";:INPUT n1
PRINT "old bmin =";bmin;" new  =";:INPUT bmin
 GOTO 12
 100 PRINT
 
 PRINT"                   d theta = ";dt
 PRINT"                  ka =";ka2
 PRINT"                theta max = ";tmax
PRINT"                n max = ";nmax
PRINT"                nmin = ";n1
PRINT"                bmin = ";bmin
GOTO 12

 
200 'compute
'compute equal log steps, 10 per decade.

FOR m = 0 TO 14
    bm(m) = EXP(-m*da)
    PRINT m,bm(m)
NEXT m
' do reduced calculation (ka)^2 w(b)
 D(0) = 1
 p(0) = 1
 wr(0) = 1.9995
 b(0) = 1

FOR n = 1 TO nmax-1
    phi = n*dt*pi/180
    x = ka2*SIN(phi)
    IF x<>0 THEN
    GOSUB 1000
   Dn = 2*J1/x
   END IF
   z = ABS(4*COS(phi)*Dn*(J0 - Dn))
 IF z>0 THEN w = (ka2*SIN(phi))^2/z
 wr(n) = ABS(w)
 b(n) = Dn^2
 NEXT n

FOR j = 1 TO nmax
    IF b(j)> bmin THEN nb = j
 NEXT j
PRINT" nb = ";nb
INPUT q$

m = 1
wdb(0) = wr(0)
FOR n = 1 TO nb
       phi = n*dt*pi/180
   IF bm(m)>b(n+1) AND b(n) >bm(m) THEN 
       wdb(m) = wr(n)
       phi(m) = phi
      m = m +1
  END IF
NEXT n

INPUT q$
mmax = m
FOR m = 0 TO mmax
   PRINT m,bm(m),wdb(m)
NEXT m
INPUT q$

'end of directivity computation	
FOR n = 0 TO nb STEP INT(nb/15)
 PRINT n,b(n),wr(n),lgcv*LOG(wr(n))
NEXT n
INPUT q$
GOTO 12


 1000   'subroutine for  J0(x) and J1(x)
  IF ABS(x)  <  3 THEN
     q=x/3:q1=q
     q2=q^2:q3=q^3:q4=q^4:q5=q^5:q6=q^6:q8=q^8:q10=q^10:q12=q^12
     J0 = 1+jc2*q2+jc4*q4+jc6*q6+jc8*q8+jc10*q10+jc12*q12
     J1 = x*(.5+ja2*q2+ja4*q4+ja6*q6+ja8*q8+ja10*q10+ja12*q12)
 ELSE
     q = 3/x
     sxq  = SQR(x)
     q1=q:q2=q^2:q3=q^3:q4=q^4:q5=q^5:q6=q^6
    fj0= jd0+jd1*q1+jd2*q2+jd3*q3+jd4*q4+jd5*q5+jd6*q6
    qj0 = x+je0+je1*q1+je2*q2+je3*q3+je4*q4+je5*q5+je6*q6
     qa1= x + jg0 + jg1*q1+jg2*q2 +jg3*q3+jg4*q4+jg5*q5+jg6*q6
     fj1 = jb0 +jb1*q1+ jb2*q2 + jb3*q3 + jb4*q4 + jb5*q5 +jb6*q6
     J0 =fj0*COS(qj0)/sxq
     J1 = fj1 * COS(qa1) / sxq
 END IF

 RETURN
 
5000  REM  ------- WRITE FILE---


 
PRINT " name THE FILE";:INPUT NN$ 

OPEN NN$ FOR OUTPUT AS #2
WRITE #2,nmax
WRITE #2,ka2
  FOR m = 0 TO mmax
    WRITE #2,bm(m), wdb(m)
  NEXT m
  
   FOR m = 0 TO mmax-1
   r = 100*TAN(phi(m+1))
   r1 = 100*TAN(phi(m))
    WRITE #2,phi(m),INT(r),INT(pi*(r^2-r1^2))
  NEXT m
 
CLOSE #2
GOTO 12
 
9000 END
 

 

